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Preface 


This document describes work performed during the past year under Grant NAG 
2-373. The work is directed towards the development of efficient computational 
methods for the solution of both inviscid and viscous transonic flow problems, with 
emphasis on problems of complex, three-dimensional geometry. Solution-adaptive 
techniques are being developed, using the multigrid method to achieve computational 
efficiency, and the work is focussed upon the aerodynamic problems associated with 
complex wing-fuselage shapes, including engine inlets. The Principal Investigator for 
the project is David A. Caughey, Professor in the Sibley School of Mechanical and 
Aerospace Engineering at Cornell University. The Grant Technical Monitor is 
Dr. W. J. Chyu 
Applied Aerodynamics Branch 
Mail Stop 227-6 
NASA Ames Research Center 

mm <*;«»■ *..94 

Further information can be obtained from: 

Professor David A. Caughey 
218 Upson Hall 
Cornell University 
Ithaca, New York 14853 
Tel: (607) - 255-3372 
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I. Introduction 

This report describes progress in research directed towards the efficient solution 
of the inviscid Euler and Reynolds-averaged Navier-Stokes equations for transonic 
flows through engine inlets, and past complete aircraft configurations, with emphasis 
on the flowfields in the vicinity of engine inlets. The research focusses upon the 
development of solution-adaptive grid procedures for these problems, and the 
development of multi-grid algorithms in conjunction with both, implicit and explicit 
time-stepping schemes for the solution of three-dimensional problems. The work 
includes further development of mesh systems suitable for inlet and wing-fuselage- 
inlet geometries using a variational approach. One portion of the research forms the 
basis of the Ph.D. thesis of Mr. Dun C. Liu, and another portion forms the basis of 
the M.S. thesis work of Mr. Ravi Iyer. Mr. Liu is supported independently by a Fel- 
lowship, and Mr. Iyer will be supported next year as a Graduate Research Assistant if 
the grant is renewed. 

II. Review of Work in Past Year 

Work in the past year has concentrated upon two-dimensional problems, and has 
been in two general areas: 

(1) the development of solution-adaptive procedures to cluster the grid cells in 
regions of high (truncation) error, and 

(2) the development of a multigrid scheme for solution of the two-dimensional 
Euler equations using a diagonalized ADI smoothing algorithm. 

The work in each of these areas will now be described briefly. 

A. Adaptive Grid Generation 

This section contains a review of the formulation of the Euler-Lagrangian equa- 
tions for grid generation, and some test results demonstrating their applicability. 

The solution of a partial differential equation can be greatly simplified by a prop- 
erly constructed grid. It is also true that a grid which is not well suited to the problem 
may lead to an unsatisfactory result, or even to instability and a lack of convergence. 
Therefore, one of the central problems in computing numerical solutions to partial 
differential equations is that of grid generation. Several advantages are obtained when 
the partial differential equation, approximated on an unequally spaced physical domain, 
is transformed to a uniformly spaced computational domain. For example, a curved 
body surface can be treated as a boundary in the computational domain, which fact 
permits easy application of the boundary conditions. The transformation methods fall 
into three basic classes: complex variable methods, algebraic methods, and partial 
differential equation methods. A concise review of these methods is given by Thomp- 
son. 1 

The partial differential equation method has been developed over the past decade. 
Thompson et al. 2 have developed a particular way of using elliptic PDE’s to generate 
two dimensional computational grids about arbitrary bodies. Brackbill and Saltzman 3 
have derived the elliptic PDE’s by minimizing an integral which characterizes undesir- 
able properties of the grid. 
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In order to determine the best grid point distribution, an a priori knowledge of 
the solution of the physical problem is desired. Unfortunately, this knowledge is una- 
vailable and only the general features of the solution may initially be understood. The 
concept of a solution adaptive grid is needed to incorporate information from the solu- 
tion in locating the grid points. With the use of an adaptive grid, the physical behavior 
of events on both large and small scales can be adequately accommodated by a vari- 
able mesh size. The effect of this is to reduce the discretization error of the numerical 
scheme by clustering the grid, in adaptation to the solution, in high gradient regions. 

For conceptual simplicity, an ability to handle problems with a number of length 
scales, and extendibility to three-dimensional problems, the method of Brackbill and 
Saltzman has been chosen for the present study. In this section, the formulation of 
Euler-Lagrangian equations for grid generation are reviewed, and some preliminary 
test results are presented to demonstrate their applicability. 

Adaptive Grid Schemes 

Brackbill and Saltzman , 3 and Brackbill 4 have extended Winslow’s 5 method for 
generating a computational mesh to include grid adaptation. In this approach a varia- 
tional technique is used to minimize a linear combination of measures of grid smooth- 
ness, orthogonality, and adaptation. The smoothness is measured by integrating the 
change of the computational coordinates over the physical domain, and can be written 
as the following form for a two dimensional problem 

h = f ((V£) 2 +(Vti) 2 )<£4 (1.1) 

D 

where dA represents a differential area in the physical space, and D denotes the physi- 
cal domain. In addition to smoothness requirement, control of mesh skewness is 
obtained if a measure of grid non-orthogonality, expressed as 

I 0 = f (V r\) 2 J 3 dA (1.2) 

D 

is minimized over the domain, where / is the Jacobian. The grid adaptation is pro- 
vided by minimizing a weighted average of the area variation over the field; an 
appropriate integral is 

I w = f wJdA (1.3) 

D 

where w is the weight function which produces the grid adaptation. If w = w(x,y) 
the solution for minimized l w corresponds to wj 2 = const, while for w = w(!;,r|), it 
corresponds to wJ = const. . 

In order to incorporate smoothness, orthogonality, and adaptation in the grid 
generator, a linear combination of the functionals is considered, 

I = I, + + foh (1.4) 

and the corresponding Euler equations are 

b i x \\ + b 2 x \t\ + b-pxm + art# + a&E , n + m fw<*l 0-5) 

and 

a l X %% + a 2*5Ti + a 3*nri + C 1 y%% + c z + c 3>Vi = fw G 2 (1.6) 
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where 

— **si "b f w i "h fa &oi 

^si f w^wi h f o^oi 

Q "h y>v Cyvi f o t~oi 

for i — 1, 2, 3. Here a si , b si , c si are the coefficients arising from minimization of the 
smoothness integral, a oi , b oi , c oi are the coefficients arising from minimization of the 
orthogonality integral, and a wi , b wi , c wi are the coefficients arising from minimization 
of the weighting-function integral. Ifw = w(x,y), 

J 2 


Gi = — 

G 2 = - 


iw 


2w dx 
J 2 d 


w 


2 w dy 


or if w = w-(£,r|), 


G 1= 


G 2 = -±( Wr{XrW ^ ) 

Normalization of the Functional 

For a uniformly spaced grid on a square, the initial value of the smoothness func- 
tional is I s ~ 2n 2 , where n represents the number of cells in both the x and y direc- 

A 

tions. For the solution adaptation functional, I w ~ ~r, where w is the mean value of 

n z 


w in the computational domain. For the orthogonality functional, I c 
the same order of magnitude, f w and f 0 are scaled to 

2 « 4 


n 


fw = /* 


w 


To have 


(1.7) 


and 


f - f — 

J O J o 2 

respectively. 

For more general geometries, it is reasonable to assume 

h ~ Xf.~ x n ~ ~ y n 
where h is the step size. With m x n cells, 

I s ~ mn 
I w ~ mnwh 4 
I n — mnh 4 


( 1 . 8 ) 


therefore 
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/ 


f w 



(1.9) 


and 



( 1 . 10 ) 


Scaling and Smoothing of the Weight Function 

To have some external control of the range of the weight function, wy is scaled, 
as proposed by Brackbill and Saltzman, 3 such that after scaling 


w ijJij = 1 • 

The range of wy is determined by the lesser of the externally set value Oq and the 
ratio of the weight function 


o = min 


_ max(w t7 ) 

13 o> — r— 

mm (Wyj 


In order that the largest zone fits into the computational domain, q, the minimum 
value of w must satisfy the equation 


q = min(w iy ) = min 


1 J_ 

ZW’-o 

ij ' 


The scaled wy is calculated from 


Wy = q 


(°- IK , l 

max(wy) 


(l.H) 


Computational experience has indicated that the convergence rate is poor when 
the weight function has very large gradients. The weight function can be smoothed by 
the following equation 


-(m+i) 

y 


(l-|l)Wjj«> + 0.25 


( 1 . 12 ) 


where |l£ 0.5 to ensure stability. 


Numerical Schemes 


All derivatives in Eqs. (1.5) and (1.6) are approximated by central differences 

x % ~ 0 . 5 ( x i+ 1 j— x i- 1 j ) 
x r\ ~ 0.5(Xj y+ j— 

and 

x $l~ x ‘+ by - ^ x i,j+ x i- ij 
~ 0.25(x i+ 1 J + { + X;_ iy_ i-X i+ IJ_ I — X,._ lj+l) 

~ x i,j+l~^ x i,j +x i,j-l 

If the residuals at the m-th iteration of this approximation to Eqs. (1.5) and (1.6) are 
denoted by R and R ^ , respectively, the Jacobi method gives 
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R[ m) = 2(b[ m) + bi m '>)&xft +l) + 2(a{ m ) + a^ m ) )Ay£j +v * 

R± m) = 2(ai m) +a^)Axtf + » + 2(c 1 <'»> + c± m) )Ay£j +1 '> 

One can solve for Ax£j +l) and A y£j +1) using Eqs. (1.13) and (1.14). The iteration is 
terminated when R f m) and R | m) are everywhere smaller than a prescribed tolerance. 
Numerical Examples 

Two cases are studied. The first deals with a model singular perturbation prob- 
lem and the second deals with the nonlifting transonic flow over an NACA 0012 air- 
foil. 

Case 1 - A Model Singular Perturbation Problem 
Consider the following convection-diffusion equation 

<)>, + V •(<!>«") — V -V <j) = 0 

where 


(1.13) 

(1-14) 


and P e 


and, 


u = -P e h + (T)h~(T) l r 
h + (T) = (1 + exp(P e (F— l))) -1 
h~{T)= (1+ exp^l-r)))" 1 

“o^o^o- If the boundary conditions are such that a steady state is reached, 

u <t> = V <() (1.15) 


<t)(r) = exp (~P e T h~(7) + log(/i + (7)) - log(/t + (0))) 

It is noted that <j> has value 1 within T= 1 and sharply drops to zero outside T= 1. 


If the weight function is chosen as 


w = 


<t> dx <j> dy 


2 

) 


with (1.15) 


w = ( u-l ) 2 + (u-j ) 2 


In order to maintain symmetry, a simple point-Jacobi method is used to solve 
Eqs. (1.13) and (1.14). For the weight function displayed in Figure A.l, the grids for 
P e = 12, r 0 = 0.25 and f w = 5 are shown in Figures A.2 and A. 3. The mesh in Fig- 
ure A.2 was computed with a simple Dirichlet boundary condition to enforce uniform 
spacing on the boundaries. That in Figure A.3 is obtained by including an orthogonal- 
ity constraint (i.e, f a = 1), and orthogonality is also enforced at the boundaries. 

Case 2 - Transonic Flow over an NACA 0012 Airfoil 

For the calculation of compressible flows with shock waves, the pressure, 3 Mach 
number, 6 or internal energy 6 have been proposed to serve as weight function for solu- 
tion adaptation. However, when using multigrid a more natural weight function has 
been proposed by Berger and Jameson, 7 which is naturally related to the local trunca- 
tion error of the scheme. This procedure has also been adopted for use in the present 
work. 




Figure A.l Weight function for circularly-symmetric test case. 



IM&! r« 


OF POOR QUALITY 



Fig. 2 


Grid 48 x 48 Iter= 500 


Figure A. 2 Mesh computed with Dirichlet boundary conditions 
and no orthogonality control. 



Fig. 3 


Grid 48 x 48 Iter= 500 


Figure A. 3 Mesh generated with orthogonality constraint. 
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Let the equation we wish to solve be represented as 

L(u) = 0, 

where L is a differential operator and u is the exact solution. Let the numerical 
approximation to this equation be represented by 

N h (U) = 0 

where N h is the numerical operator on the finest grid with step size h, and U is the 
numerical approximation to u (i.e., is the exact solution to the difference equations'). 
Since 


the local truncation error is 


u = U + AU , 


N h (u ) - h n D , 


where D contains derivatives of u, and n = 2 for a method which is second order 
accurate. The goal of the adaptive procedure is to find where D is large and to reduce 
the mesh spacing h there. To estimate this error we use the relation on a coarser grid 
having twice the mesh spacing of the fine grid 

N 2h (u) - (2 h) n D . 


If AC/ is small, 


N Zh (U) - (2 h) n D , 


and interpolating back to the fine grid gives 

/2*(tf”(ff)) ~ h*D . 
For steady state calculations, this relation is modified as 


h n D - I h 2h (N 2h (U))-N h (U) . (1.16) 

If t ^ ie calculation is started with a uniform mesh distribution, regions where 
l^(^ 2A (C/))-iV A (C/)|is large indicates D is large and h needs to be reduced. 

Since the solution minimizing I w for w = w(^,T[) corresponds to wJ = const . , the 
weight function is taken as 


|/2A(^ 2A (C/))-iV /, (C/)| 

J 


(1.17) 


with the adaptation process terminating when h n D reaches a constant. 

The 64x 16 O-grid used in Flo52 in calculating flow field over an NACA 0012 air- 
foil at Mach number 0.80 and zero angle of attack is taken as a test case for applying 
the solution adaptive algorithm. The residual calculation, based on the root mean 
square of the term dpldt after 50 multigrid cycles, indicates that the regions of largest 
truncation error are located near the leading and trailing edges of the airfoil. With the 
weight function defined as in Eq. (1.17), the Jacob i method is used to solve Eqs. 
(1.13) and (1.14), with the far field boundary points held fixed and the inner boun- 
dary points allowed to move along the profile surface to satisfy grid orthogonality. 
After 100 iterations the new grid is shown in Figures A.4 - A.7. The pressure 



Old Mesh from flo52 
Grid 64 x 16 

Window Size XI = -0.0 X2 = 200.0 


Fig. 4 New Mesh after Adaption 
Grid 64 x 16 Iter= 100 
Window Size XI = -0.0 X2 = 200.0 


Figure A. 4 Original and adaptive meshes for transonic airfoil problem. 




Old Mesh from flo52 
Grid 64 x 16 

Window Size XI = 97.0 X2 = 103.0 



Fig. 5 New Mesh after Adaption 
Grid 64 x 16 Iter= 100 
Window Size XI = 97.0 X2= 103.0 


Figure A. 5 Original and adaptive meshes for transonic airfoil problem. 




Old Mesh from flo52 
Grid 64 x 16 

Window Size XI = 97.8 X2 = 98.4 



Fig. 6 New Mesh after Adaption 
Grid 64 x 16 Iter= 100 
Window Size XI = 97.8 X2 = 98.4 


Figure A. 6 Original and adaptive meshes for transonic airfoil problem. 



Old Mesh from flo52 
Grid 64 x 16 

Window Size XI = 101.6 X2 = 102.2 


Fig. 7 New Mesh after Adaption 
Grid 64 x 16 Iter= 100 
Window Size XI = 101.6 X2 = 102.2 



Figure A. 7 Original and adaptive meshes for transonic airfoil problem. 
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distribution and convergence history after the calculation is continued for another 50 
multigrid cycles on this new grid are shown in Figures A. 8 and A.9, respectively. It is 
worth noting that after the same number of multigrid cycles the solution adaptive grid 
has reduced the largest truncation errors at the leading and trailing edges, and the 
standard deviation of the truncation errors on the original grid, by almost fifty per 
cent. 

B. Diagonal Implicit Alternating Direction Multigrid Scheme 

In order to circumvent the problem of slow convergence of explicit methods on 
grids containing cells of high aspect ratio, it seems attractive to develop implicit alter- 
nating direction multigrid methods. These methods have also been studied recently 
by Jameson and Yoon. 8 In order for the implicit method to be an effective smoothing 
algorithm when used in conjunction with the multigrid algorithm, it is important to 
include an accurate representation of the dissipative terms. Since these usually include 
fourth-differences to maintain high accuracy, their inclusion in the implicit operator 
requires the inversion of block pentadiagonal systems for each one-dimensional factor. 
To avoid the high cost of solving block pentadiagonal systems, the equations in the 
present scheme are first diagonalized at each point using a local similarity transforma- 
tion, following Chaussee and Pulliam.^ This has the effect of decoupling the equations, 
and requiring the solution of four (in two-dimensional problems) scalar pentadiagonal 
equations for each factor. The resulting method has good high-wavenumber damping, 
so is a good smoothing algorithm for use in conjunction with the multigrid method. It 
is also computationally efficient because of the need to solve only scalar systems. 
Additional computation is required to calculate the local similarity transformations, 
and to perform matrix multiplies of the residual, and of the intermediate and final 
corrections, but this is a small fraction of the work required to solve the block sys- 
tems, and also can be vectorized to further reduce the required CPU time. 

The Euler equations of inviscid compressible flow can be written (in two space 
dimensions) as 

drt/dt + df^/dx + d]?/dy — 0 , (2.1) 

where 

# = (P, P«, pv, e} T , (2.2a) 

is the vector of dependent variables, and 

/= {pu, p u 2 +p, p uv, (e+p)u} T , (2.2b) 

{ pv, p«v, pv 2 +p, (e+p)vf, (2.2c) 

are the flux vectors in the x and y coordinate directions, respectively. Here, p and p 
are the fluid density and pressure, «,v are the velocity components in the x and y 
directions, and e is the total energy per unit volume. The pressure is related to the 
total energy per unit volume e by the equation of state 

P - (Y-l){e - P(k 2 + v 2 )/2}, (2.3) 

where y is the ratio of specific heats. 





0. 50. 100. 150. 200. 250. 300. 


NACA0012 

Mach 0.800 Alpha 0. 

Residl 0.102e+01 Resid2 0.639e-04 

Work 131.48 Rate 0.9290 

Grid 64x 1 6 



0. 50. 100. 150. 200. 250. 300. 


Fig. 9 

Mach 0.800 Alpha 0. 

Residl 0.102e+01 Resid2 0.849e-05 

Work 132.81 Rate 0.9157 

Grid 64x16 


Figure A. 9 Convergence histories for original and adaptive grid solutions. 
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The algorithm has been implemented within a finite-volume framework to allow 
the treatment of essentially arbitrary geometries. Under an arbitrary nonsingular 
transformation of dependent variables the equations can be written 

9^ /dt + dft/dZ, + dCr/dT[ = 0 , (2.4) 

where W = hrt is the transformed dependent variable and 

f= {phU, phUu+y^p, phUv-x^j, (e+p)hU} T , (2.5a) 

S - {p hV , phVu-yy>, p hVv+xy>, ( e+p)hV } T , (2.5b) 

are the transformed flux vectors. Here h = xy^-xy^ is the determinant of the Jaco- 
bian of the transformation (which corresponds to the cells area), and U and V are the 
contravariant components of the velocity given by 


h{U,V} T = 


y n “ - 

— y^u + x^v 


( 2 . 6 ) 


In a finite-volume method, the spatial derivatives are represented by approximating 
the net flux across the faces of each mesh cell, using constant values of the velocities 
to evaluate the flux across each face. In the present method, the dependent variables 
are defined at the cell centers; alternatively, these values can be considered as aver- 
ages for the cell. The value on each face is taken to be the average of the cells sharing 
the face, so that for example 

(hU) i+i/it j = V2{(yij + i A — yij_ J/4 ) ( u i+ + U;j) - 

x iJ-ti(v i+hj + V;j)} (2.7) 

This approximation is equivalent to a centered difference scheme which is 
second-order accurate in the mesh spacing in the physical domain when the mesh is 
smooth. In order to prevent the decoupling of the solution at odd and even num- 
bered cells in the grid, dissipative terms must be added. Following Jameson, 10 the 
dissipative terms are constructed as an adaptive blend of second- and fourth- 
differences. As pointed out by Jameson, the fourth-difference terms are necessary if 
the scheme is to converge to a steady state, while the second-difference terms are 
necessary to prevent excessive oscillation of the solution in the vicinity of shock 
waves. 

The difference approximation, including the dissipative terms, can be written 

d(Wij)/dt + Qrfij - Drfij = 0 , ( 2 . 8 ) 

where Q is the operator corresponding to the differences introduced by approxima- 
tions of the form of Eq.(2.7) for the fluxes, and D is an operator corresponding to the 
dissipative terms, which will now be defined. The dissipative operator can be written 

£>v? = D y? + . (2.9) 


where 


= S^dyit), 

D T\rf= V ^), 


(2.10a) 

(2.10b) 
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where 8^ and S n are central difference operators spanning a single mesh cell, and 

{e (2) 5^ - e (4) 5|}>^ . (2 .H) 

The coefficients e (2) and e (4) are adapted to the solution and are defined as follows. 
Let 


then 


and 


v /,y - I (Pi+l,j ~ 2Pi,j+Pi - 1 j)/ (Pi+ l y j + 2pij+p i _ l j ) | , 


e ‘+kj = * (2) max(v /+1J , Vij) 


4\j = max(0, hlJr ' - A ; J k (4) 
At 


e/+ky) 


( 2 . 12 ) 


(2.13a) 


(2.13b) 


where At* is the time-step corresponding to unit Courant number for the cell, and k ( 2) 
and k (4) are constants (chosen to be 1 and 1/64, respectively, in the present study. 
This scaling of the dissipative terms with 1/Af* makes them proportional to the propa- 
gation speed of that characteristic which limits the time step. The dissipative terms 
corresponding to the operator D ^ are similarly defined. This tailoring of the dissipa- 
tive coefficients simultaneously activates the second-difference terms and turns off the 
fourth -difference terms near shock waves, allowing them to be captured with little, or 
no, overshoot. 

In smooth regions of the flow, the equations corresponding to Eqs.(2.4) plus the 
dissipative terms just described can be written in the quasilinear form 

dW' /dt + A 3W^ Id^A BdW Id TJ = 

a{(e (2) a>wa£- e (4) a 3 w/3^ 3 )}/a| + 

3{(e (2) 3vwari - e (4) a 3 w/aTi 3 )}/ar| , (2.14) 

where A = {dA/dtf} and B = {dCr/aw^} are the Jacobians of the flux vectors with 
respect to the solution. These matrices are given by Chaussee and Pulliam. 9 

The basis of Alternating Direction Implicit (ADI) methods is to approximate the 
spatial derivatives as weighted averages of differences taken at the old and new time 
levels, linearizing the changes in the flux vectors using the Taylor series expansions in 
time 

ft,j l = ?i.j + KA^li + 0(At 2 ), (2.15a) 

^ij 1 = + B'jAtfTj + O {At 2 ), (2.15b) 

where 

AV?. n . = U^t 1 _ 

is the correction to be added to the solution in the i,j cell in going from the n to the 
n+ 1 time level. This gives a scheme of the form 

{/ + BAt [Aij% + Bij5 n - e/ 2 >(5| + S 2 ) + e/ 4 >(8( + 5 4 )]f Atf? tj = 

- Atfyfij + 8 n g itj + eg>(Sf + 5 2 )tf - eW{8$ + 5 4 )v?]}", (2.16) 
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where 0 is a parameter determining the degree of implicitness of the scheme (with 
0=1 corresponding to a fully implicit scheme. 

Note that for simplicity, Eqs.(2.16) are written as if the coefficients of the and rt- 
difference dissipation terms are the same, but they are, of course, different, depending 
as they do on the second derivative of the pressure in the appropriate coordinate direc- 
tion. 


Finally, to make the scheme computationally efficient, the operator on the left- 
hand side of Eq.(2.16) is approximated as the product of two one-dimensional factors 
to give 

{/ + 0Af [A" 7 -5| - e/,?5f+ e/j>8£|} 

{/ + 0Ar[^ 7 5 n - effS* + e^5 n 4 ]}A^ = 

- Atfyfij + 8 n & id + e/,5>(8| + 8 2 )v? - e/j>(8^ + 8 n V]} n - (2.17) 


The solution of Eqs.(2.17) is achieved in two steps. First, the intermediate correction 


AV/’j ={1 + 64 


(2.18a) 


is determined by solving 

{/ + 0A tlAWt ~ tQH + e^8(]}A^ y = itij , (2.18b) 

where R t j is the residual vector corresponding to the right hand side of Eqs.(2.17). 
Thus, to advance the solution one time-step, Eqs.(2.17) require the inversion of one 
block pentadiagonal system along each line of constant r|, followed by the inversion of 
one block pentadiagonal system along each line of constant The size of the blocks 
is {Ax A) for this two dimensional problem; for a three-dimensional problem, there are 
three factors to be inverted, and the blocks are (5x5). 


The scheme described by Eqs.(2.17) would be reasonably efficient if it were not 
necessary to add numerical dissipation to stabilize the central-difference approximation. 
Even so, if only second-difference terms were added (i.e., if e (4) =0), it would still be 
necessary only to solve block tridiagonal systems. As pointed out above, however, the 
inclusion of fourth-differences is essential if the solution is ultimately to converge to a 
steady state, and it is important to treat them implicitly if the solution is to converge 
rapidly. This can be done in a straightforward manner, following the development 
above, but leads to a requirement to solve block pentadiagonal systems for each factor. 
This requires approximately twice the computational labor of the block tridiagonal 
inversion, and begins to become computationally prohibitive. 

An alternative is to diagonalize the equations at each mesh point, yielding a 
decoupled set of equations which can be solved using a scalar pentadiagonal solver. 
This requires approximately one-quarter the computational labor of the block pentadi- 
agonal solution (and requires, in fact, only about half the computation required to 
solve the block tridiagonal systems). Let M A and M B be the modal matrices of the 
Jacobians A and B , so that Q a 1 AQ a = A a and Qb 1 BQ b = A B are diagonal matrices 
whose elements are the eigenvalues of A and B , respectively. At each point the equa- 
tions can then be written 


{/ + 0At[AjJS ? - e (2) 8| + e (4) S£j} 
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*> 


Qa 1 Qb {/ + ©At [A fl B S n - e (2) 5^ + e< 4) 8 4 ]}Av£,- = 

- Ar QZ'fyftj + 8 + eg>(8f + 8 n V - e<?(S 5 4 + 8j)0]r (2.19) 

where 

A#/|; = Q b AV>?j , 

i.e., AV? j is the change in the characteristic variable corresponding to the Jacobian of 
the t| flux vector. The solution procedure is similar in structure to that of the block 
scheme. First, the intermediate correction 

A V7j = Qa 1 Qb {/ + 9*ttAS\ - e (2) 5 2 + £ (4) 8 4 ]}AV^j (2.20a) 

is determined by solving the equations 

{/ + 8Ar[A£S 4 - e (2) S|+ e (4) 8|]}Av7j = Qa'^j , (2.20b) 

then the correction itself is determined by solving Eqs.( 2.20a). As described above, 
the solution of Eqs.(2.19) requires the determination of the modal matrices of the A 
and B Jacobians (and their inverses) at each point, a matrix multiply of the residual 
vector, the solution of four scalar pentadiagonal systems along each line, another 
matrix multiply of the intermediate correction vector, a second set of four scalar pen- 
tadiagonal solutions along each ri-line, and a final matrix multiply to extract the 
corrections to the primitive variables. The cost of computing the elements of the 
modal matrices is about the same as that of computing the elements of the Jacobian 
matrices in the block method, and even with the extra matrix-vector multiplies, the 
diagonalized procedure is still considerably more efficient than solution of the block 
pentadiagonal systems. The advantage of the diagonal scheme would be even greater 
on a vector machine, since the additional matrix-vector multiplies are vectorizable. 

The incorporation of the scheme within the multigrid algorithm is straightfor- 
ward, following the procedure developed by Jameson. 11 The procedure will not be 
further described here. Since Jameson found it desirable to use only a fixed- 
coefficient, second-difference form of the dissipation on coarser grids, additional 
efficiency with the present formulation can be achieved by using scalar tridiagonal 
solvers on the coarser grid levels. 

The algorithm described above has been coded for the problem The diagonal 
implicit scheme has been coded for the problem of transonic flow past an airfoil, 
including the blended second- and fourth-difference dissipation terms and the mul- 
tigrid implementation. Results have been obtained for a variety of airfoils and frees- 
tream conditions, to verify the accuracy of the basic algorithm. Here, the results of a 
single computation will be presented to verify the efficiency of the multigrid imple- 
mentation. 

The calculation was performed on a C-grid containing 192 x 32 grid cells in the 
wraparound and normal directions, respectively. The farfield boundary of the grid is 
located approximately 50 chords upstream and downstream, and 100 chords laterally, 
from the airfoil, and the most elongated cells have aspect ratios of approximately 168 
and 114 in the r\ and S, directions, respectively. The grid structure in the vicinity of 
the airfoil surface is shown for a NACA 0012 profile in Figure B.l. The calculations 
were performed on an FPS AP-264 Scientific Processor attached to an IBM 3084 host. 




NACA 0012 Airfoil 
Grid 192 x 32 


Figure B.l Computational grid in vicinity of airfoil surface. 
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A typical calculation (consisting of 100 work units) requires approximately 3.0 minutes 
for the decoupled pentadiagonal implicit scheme, or 8.4 minutes for the tridiagonal 
block implicit scheme, demonstrating the relative computational efficiency of the diag- 
onalized scheme. 

The airfoil surface pressure distribution for the case of the flow past the NACA 
0012 airfoil at a freestream Mach number of 0.80 and 1.5 degrees angle of attack is 
presented in Figure B.2, showing a strong shock on the upper surface, and a weak 
shock on the lower surface of the profile. 

Convergence results are presented in Figure B.3 for this case. The solid lines in 
Figure B.3 represent the convergence history for the present diagonal implicit mul- 
tigrid scheme for this transonic, lifting case, while the broken line represents the con- 
vergence history for a single-grid implementation of the diagonal implicit algorithm. 
The logarithm of the average residual (average over all the cells of the absolute value 
of Ap / At) and the total number of cells in which the local Mach number is supersonic 
are plotted vs. work units. The latter quantity is a measure of the global convergence 
of the solution. The Courant number for the implicit multigrid calculation is 8., while 
that for the single-grid calculation is 12. For both calculations, local time-stepping is 
used. (That is, the time step in each cell is adjusted to correspond to a fixed Courant 
number). The implicit scheme remains stable at higher Courant numbers, but the 
overall convergence rate is degraded because of the poorer high-wavenumber damping 
at the larger values of time step. The average residual in the multigrid calculation has 
been reduced by almost seven orders of magnitude in only 100 multigrid cycles 
(corresponding to approximately 200 Work units) and the number of supersonic cells 
has converged to its final value after only 50 multigrid cycles (corresponding to about 
100 work units). The single-grid residual has been reduced by less than three orders 
of magnitude in the same 200 work units, and the number of supersonic cells is still in 
error by more than 7 per cent. 

Work is continuing on the development of the scheme, with particular attention 
being given to the construction of the dissipative terms, and the performance of the 
scheme on more highly stretched grids. In addition, the algorithm is being applied to 
compute the flow in a supersonic inlet. A computational grid suitable for the multigrid 
scheme has been developed using a modification of the procedure of Caughey and 
Jameson; 12 a sample of the grid is shown in Figure B.4. Attempts are also being made 
to modify the GRAPE code 1 ^ for this class of geometries in order to provide greater 
control over local mesh spacing. Work on adapting the diagonalized ADI multigrid 
algorithm for this problem is currently underway. 
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-1.60 h 
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NACA 0012 Airfoil 
Mach 0.800 Alpha 1.250 
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Figure B.2 Airfoil surface pressure distribution. 
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Figure B.3 Convergence histories with and without multigrid. 
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Figure B.4 Computational grid for inlet flow field. 
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